Zero-temperature transition and correlation-length exponent of the frustrated XY 

model on a honeycomb lattice 
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Phase coherence and vortex order in the fully frustrated XY model on a two-dimensional honey- 
comb lattice are studied by extensive Monte Carlo simulations using the parallel tempering method 
and finite-size scaling. No evidence is found for an equilibrium order-disorder or a spin/ vortex- glass 
transition, suggested in previous simulation works. Instead, the scaling analysis of correlations of 
phase and vortex variables in the full equilibrated system is consistent with a phase transition where 
the critical temperature vanishes and the correlation lengths diverge as a power-law with decreas- 
ing temperatures and corresponding critical exponents v p h and v v . This behavior and the near 
agreement of the critical exponents suggest a zero-temperature transition scenario where phase and 
vortex variables remain coupled on large length scales. 
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I. INTRODUCTION 



Two-dimensional frustrated XY models 1,2 , in which 
vortices form a dense incommensurate lattice, have 
attracted considerable interest as a possible two- 
dimensional vortex glass without quenched disorder 5-7 
or a structural glass of supercooled liquids^. In su- 
perconducting systems with pinning described by XY 
models, frustration can be introduced by applying an 
external magnetic field as in periodic arrays of Joseph- 
son junctions^— , superconducting wire networks^ 3 - - — 
and superconducting thin films with a periodic pattern 
of nanoholesi^. The frustration parameter / sets the av- 
erage density of vortices in the lattice of pinning centers 
and can be tuned by varying the strength of the external 
field 2 . Depending on the structure of the lattice of pin- 
ning centers and the value of /, a commensurate vortex 
lattice is favored in the ground state, which leads to dis- 
crete symmetries in addition to the continuous symmetry 
of the phase variables characterizing the superconducting 
order parameter. In this case, the phase transitions and 
resistive behavior of the system are reasonably well un- 
derstood for simple low-order commensurate phases such 
as / = 1/2, on a square latticed, and / = 1/3, on a tri- 
angular latticed of pinning centers. However, when the 
vortex lattice is incommensurate with the pinning cen- 
ters, both the nature of the equilibrium phase transition 
and of the low-temperature state in the thermodynamic 
limit is less clear. This is the case of Josephson-junction 
array s 10 ' 12 and superconducting wire network s 13 ' 14 on a 
square lattice, described by the frustrated XY model with 
irrational /, which has been extensively studied by vari- 
ous methods 3 -^"— ' 18 ' 19 . Another interesting physical re- 
alization but much less investigated so far, should occur 
in Josephson-junction arrays and superconducting wire 
networks 1 - 5 - on a honeycomb lattice, and superconducting 
films with a triangular pattern of nanoholes^, which can 
be described by the fully frustrated XY model (/ = 1/2) 
on a honeycomb lattice. 



In early Monte Carlo (MC) simulations of the fully 
frustrated XY model on a honeycomb lattice^ 7 -, a phase 
transition at a nonzero temperature in the Koterlitz- 
Thouless universality class was suggested, based on the 
saturation of the specific-heat peak with increasing sys- 
tem sizes and the apparent jump in the helicity modu- 
lus. On the other hand, calculations of the the spin-glass 
order parameter by MC simulations 5 - suggested that, in- 
stead of an order-disorder transition, a spin-glass transi- 
tion should take place approximately at the same tem- 
perature. In contrast, from MC simulations of a simi- 
lar model in the vortex representation 2 ^,, it was argued 
that only a crossover region rather than an equilibrium 
phase transition should occur at any nonzero tempera- 
ture, since the energy of domain-wall excitations which 
disorder the ground state was found to approach a finite 
constant for increasing system sizes. However, the free 
energy of domain-wall excitations in the frustrated XY 
model obtained analytically in the phase representation 21 
was found to increase linearly with the system length, 
but with an extremely small numerical coefficient. As 
a consequence, although vortex ordering could be possi- 
ble in the thermodynamic limit, domain-wall excitations 
would only be negligible for system sizes which are much 
beyond the ones studied numerically or even experimen- 
tally. It is unclear if the behavior observed in the earlier 
numerical work s 5 ' 17 could be a signature of such vortex 
ordering or the effect of slow dynamics which prevents to 
observe the true equilibrium behavior. Given these con- 
flicting results, it should be of interest to further investi- 
gate the fully frustrated XY model using a MC method 
which can insure full equilibration and a scaling analysis 
of both phase and vortex correlations. 

In this work we study phase coherence and vortex or- 
der in the fully frustrated XY model on a honeycomb lat- 
tice, using the parallel-tempering method (exchange MC 
method) 22 ' 23 to obtain equilibrium configurations of the 
system. This method has been shown to reduce signifi- 
cantly the long equilibration times in glassy system s 22 ' 24 . 
To study the equilibrium phase transitions we use numer- 
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ical data in the temperature regime in which full equili- 
bration can be insured and employ a scaling analysis of 
the correlation lengths of the phase and vortex variables 
to extrapolate to the low-temperature and large-system 
limits. No evidence of an equilibrium order-disorder 
phase transition or even a spin/ vortex-glass transition is 
found at nonzero temperatures. Our results, however, are 
consistent with an equilibrium zero-temperature transi- 
tion, where the critical temperature vanishes (T c = 0) 
and the correlation lengths diverge as a power-law with 
decreasing temperatures and corresponding critical expo- 
nents v p h and v v . This transition has important conse- 
quences for the resistivity behavior and current-voltage 
scalin g 6 ' 31 ! 32 at low temperatures in the superconducting 
systems described by the frustrated XY model. More- 
over, the near agreement of the critical exponents esti- 
mated numerically suggests a T c = transition scenario, 
where phase and vortex variables remain coupled on large 
length scales. This is in contrast with the T c — tran- 
sition in the frustrated XY model on a square lattice 
at irrational frustration, where a decoupled scenario was 
found recently^. 



II. MODEL AND MONTE CARLO 
SIMULATION 

We consider a uniformly frustrated XY model de- 
scribed by the Hamiltonian 2 

h = -jJ2 cos (^ - h - (!) 

<ij> 

where 9i is a phase variable defined at the sites i of a 
two-dimensional honeycomb lattice (Fig. 1), representing 
the local angle of the XY spin with an arbitrary fixed 
direction. J > is a uniform ferromagnetic coupling 
and Aij is a gauge-invariant quantity constrained to be 
Ylij Aij = 27t/ around each plaquette of the lattice. The 
parameter / controls the frustration of the system. For 
the fully frustrated case considered here, / = 1/2. In 
the calculations we choose a gauge where = 2nfni/2 
on the bonds along the horizontal rows numbered by the 
integer rii and A^ — on the vertical bonds of the lattice. 

As a model of an array of superconducting grains cou- 
pled by Josephson junctions, the phase Qi in Eq. (TT]) 
corresponds to the phase of the superconducting order 
parameter of the grains, J is the Josephson coupling be- 
tween grains and A^ is the line integral of the vector 
potential A due to an external magnetic field B = V x A 
applied perpendicular to the array. The magnetic flux in 
each plaquette in units of the flux quantum $ = hc/2e 
can be written as 27r/ with the frustration parameter / 
corresponding to the number of flux quantum per pla- 
quette. 

In the Monte Carlo simulations we use the parallel tem- 
pering method 22 to obtain equilibrium configurations. In 
this method, many replicas of the system with different 




FIG. 1: Honeycomb lattice and the coordinate axes. 

temperatures are simulated simultaneously and the cor- 
responding configurations are allowed to be exchanged 
with a probability satisfying detailed balance. The ex- 
change process allows the configurations of the system 
to explore the temperature space, being cooled down 
and warmed up, and the system can escape more eas- 
ily from metastable minima at low temperatures. Full 
equilibration can be insured within reasonable computer 
time in systems of sufficiently small sizes22. Without 
the replica exchange step, the method reduces to con- 
ventional MC simulations at different temperatures. We 
performed MC simulations using the heat-bath algorithm 
for each replica, simultaneously and independently, for a 
few MC passes. Periodic boundary conditions were used 
on lattices with rhombic geometry (Fig. 1) of side L, 
containing N — 21? sites. Exchange of pairs of replica 
configurations at temperatures Ti and Tj and energies Ei 
and Ej is attempted with probability mm(l, exp(— A)), 
where A = (1/1; - l/Tj)(Ej -EA, using the Metropolis 
scheme. The equilibration time to reach thermal equi- 
librium can be measured as the average number of MC 
passes required for each replica to travel over the whole 
temperature range. We used typically 4 x 10 6 MC passes 
for equilibration with up to 120 replicas and 1.6 x 10 7 
MC passes for calculations of average quantities. 

III. CORRELATION LENGTH AND SCALING 
ANALYSIS 

For the study of phase coherence, we consider the over- 
lap order parameter— of the phase variables defined as 
1ph{j) = exp(i0j — 70 2 ), where 1 and 2 denote two ther- 
mally independent copies of the system, with the same 
parameters J and / in the model of Eq. [TJ At high tem- 
peratures, where each copy is thermally disordered, the 
correlation function 

C ph {r) = — ^2 < q P h{j)q P h{j + r) > (2) 

3 

is short ranged, decaying exponentially with r, while at 
low temperatures it is long ranged if an ordered phase 
or a glassy phase exits. The latter possibility was sug- 
gested in ref. d based on MC simulations and the glass 
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phase was characterized by the behavior of the Edwards- 
Anderson order parameter qEA, which corresponds here 
to < EjQph(j)/N >. 

The correlation length in the finite system £ p h(L,T) 
can be obtained a o 24,25 



Zph(L,T) = 



1 



2sin(fc /2) S ph (k 



^ Sph(0) _ jy/2 



(3) 



where S p h(k) is the Fourier transform of C p h{r) and k 
is the smallest nonzero wave vector in the finite system. 
This expression for the correlation length £(L,T), which 
is both temperature and size dependent, can be obtained 
from the correlation length £(T) in the large system limit, 
£(T) 2 = -gm y U=cb as a finite-difference approxi- 
mation taking also into account the lattice periodicity 24 . 
Thus, £(L,T) tends to £(T) for large L, when corre- 
lations are short ranged. If there is long-range order, 
£(L,T) ~ L 1+d / 2 , while if there is only algebraic or- 
der £(L,T) ~ L. Note that the correlation length de- 
fined in Eq. in terms of the overlap order parameter 
Qph(j), may have a different magnitude from the correla- 
tion length defined in terms of a single copy of the system 
exp(i#j). However, they should have the same leading 
scaling behavior near the transition. 

Analogous expressions are used to determine the cor- 
relation length for vortex variables £ V (L,T), 



UL,T) = 



1 



S v {0) 



2sin(fc G /2) S v (k 



-1) 1/2 , 



(4) 



in terms of the overlap order parameter q v (p) = v p v p of 
the net vorticity v p = n p — f . The vorticity n p is defined 
as ?ip = — Oj — Aij) /2tt, where the summation is 

taken around the elementary plaquette located at sites p 
of the dual lattice, formed by plaquette centers, and the 
gauge-invariant phase difference is restricted to the inter- 
val [— ir,n]. For the fully frustrated case, / = 1/2, the 
vortex variables v p = ±1/2 are Ising-like variables and 
equivalent to the chirality variables originally introduced 
by Villain 1 . 

Finte-size scaling can be used to extrapolate the behav- 
ior of the system to the large-system limit and temper- 
atures near the transition. In the scaling analysis of the 
correlation length 2 ^, we consider the dimensionless ra- 
tio £(L,T)/L which, for a continuous transition, should 
satisfy the scaling form 



t{L,T)/L = G{{T-T c )L 1 l») 



(5) 



where T c is the critical temperature and v is the critical 
exponent of the power-law divergent correlation length 
f(T) oc \T - T c \- V . The scaling function G(x) has the 
properties: G(0) = C, a constant, and G(x) — > x~ v as 
x — > co. This scaling form implies that data for the scaled 
correlation length £(L,T)/L as a function of tempera- 
ture, for different system sizes L, should come together 
for decreasing temperatures and cross at the same tem- 
perature T = T c . In addition, the data should splay out 



for different system sizes with slopes determined by the 
critical exponent v. If the data satisfy the scaling form 
of Eq. §5§ then we can infer that the correlation length 
diverges as a power-law £(T) oc |T — T c \~ u and estimate 
the critical exponent v. A very large value of v would be 
an indication that the correlation length may diverge ex- 
ponentially, £(T) oc e c /' T_Tc ' " , rather than as a power 
law. In particular, for the standard (unfrustrated) XY 
model, which is critical at temperatures equal and below 
T c , curves of £(L,T)/L as function of decreasing tem- 
peratures for different system sizes merge 2 ^ near T c and 
remain L- independent for T < T c . . 



IV. RESULTS AND DISCUSSION 

We study the behavior of the correlations of phase 
and vortex variables in thermal equilibrium, which im- 
plies that we should only use the numerical data ob- 
tained from the MC simulations in the temperature range 
where full equilibration is achieved. To check equilibra- 
tion, we follow the trajectory in the temperature space 
of a replica starting at the highest temperature, where 
the system equilibrates fast even without the replica ex- 
change process^. In equilibrium, this replica should be 
able to explore the whole temperature range. For the sys- 
tem sizes studied, L > 24, equilibration was only possible 
for temperatures above Tt « 0.11, despite the improve- 
ment of the parallel tempering method. Below Tf, the 
configurations of the replicas at different temperatures 
cannot be warmed up and cooled down. Thus Tf can be 
regarded as a freezing temperature, below which the sys- 
tem remains trapped in metastable configurations within 
the available time scale of the present simulation. It is 
interesting to note that Tf agrees with the apparent glass 
transition temperature T g = 0.11 observed in earlier MC 
simulations 5 . However, as will be shown in the following, 
Tf does not correspond to the critical temperature of an 
equilibrium glass transition. Instead, it should be re- 
garded as the characteristic temperature of a dynamical 
freezing transition. 

Figures 2 and 3 show the temperature dependence of 
the scaled correlation length in the ^-direction for phase 
variables, £ p h/L, and for vortex variables, £ V /L, in the 
temperature range (T > Tf) where full equilibration was 
possible and for different system sizes. Both quantities 
increase faster on lowering the temperature as the system 
size L increases indicating a divergent length scale for de- 
creasing temperatures. However, for fixed temperature 
they decrease with L even at the lowest available tem- 
perature. Therefore, the curves for the different system 
sizes do not cross or merge at a common temperature, 
as would be expected from the scaling form of Eq. |5]), 
if a transition occurs in this temperature range. This is 
more evident in the insets of Figs. 2 and 3, where lines 
joining the data for different system sizes are presented 
as a log-log plot of L/£(L, T) versus T. 

For the sake of comparison to the behavior of the vor- 
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FIG. 2: Scaled correlation length of phase variables £ p h/L for 
different system sizes L. Inset: Log-log plot of L/£ p h versus 
T for the corresponding system sizes L. 



FIG. 4: Scaled correlation length of vortex variables £ V /L 
for different system sizes L, when the frustration parameter 
is / = 1/3. 
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FIG. 3: Scaled correlation length of vortex variables £ V /L for 
different system sizes L. Inset: Log-log plot of L/£ v versus T 
for the corresponding system sizes L. 



tex correlation for the fully frustrated case (/ = 1/2) in 
Figs. 3, we show the results of additional calculations 
when the frustration parameter is / — 1/3, in Fig. 4. 
In this hexagonal vortex pattern commensurate 

with the honeycomb dual lattice is possible in the ground 
state and an order-disorder phase transition is expected 
at finite temperature from symmetry arguments. From 
the earlier worki 7 -, this transition takes place at a criti- 
cal temperature T c ~ 0.23. Indeed, the curves in Fig. 4 
cross for different system sizes at a common temperature 
T c = 0.226(1), which is compatible with this estimate. 

For the fully frustrated case, / = 1/2, the lack of in- 
tersection of the curves of £, p h/L and £«/£ for different 
system sizes at a common temperature in Figs. 2 and 3, 
suggests that phase coherence and vortex order, or even 
glass-like order, can only occur at T << Tf, which is not 
accessible in our calculations, or else only at T = 0. The 
latter scenario corresponds to a phase transition where 
T c = and the correlation length £(T) is finite at any 
nonzero T but diverges when approaching T = 0. To ver- 



ify the possibility of such zero-temperature transition, we 
first consider the scaling of the total correlation function 

X defined as x P h = J2 r C ph( r ) and Xv = J2 r C v( r ) for 
phase and vortex variables, respectively. In the absence 
of finite-size effects, \ should diverge a s 26 ' 27 

where A and B are constants and 7 is a critical exponent. 
B represents possible background corrections to scaling. 
Using data for a large system size, where finite-size ef- 
fects are negligible, we can then determine T c from the 
best numerical fit according to this scaling form. Using 
a least-square we have obtained the error estimate 
cr(T c ) when fitting Log(x — B) against Log(T — T c ) to a 
straight line, assuming different values of T c . The error 
estimate is defined as 

N p 
^ i— 1 

where N p is the number of data points {xi, yi} and f(xi) 
are the values of the fitting function at the correspond- 
ing data points. Figs. 5 and 6 show the dependence of 
the error a on T c , for phase and vortex variables, respec- 
tively, for L = 40. In both cases, T c — gives the lowest 
error, which suggests that a zero-temperature transition 
is possible. The insets in the Figures show the corre- 
sponding best fits, which provide estimates of the critical 
exponents ^ p h = 1.57(1) and -f v — 1.14(1). 

The T c = scenario can be further verified from 
the finite-size scaling analysis of the correlation length 
£(L,T). In this case, the data for £(L,T) should sat- 
isfy the finite-size scaling form of Eq. <j5j) with T c = 0. 
The best data collapse is obtained by adjusting a single 
parameter, providing an estimate of this critical expo- 
nent. We use the following procedure, which is a simpli- 
fied version of more general methods of measuring data 
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FIG. 5: Variation of the fitting error a against T c , for the 
total phase correlation Xph according to Eq. [5] with B = 2.6. 
Inset: log-log fit assuming T c = 0. 
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FIG. 7: Scaling plot according to Eq. ([5]), assuming T c = 
0, for the phase correlation length £ p h, giving the estimate 
v P h = 1.29(15). Dashed line is a fit to the data with Eq. (JHJ) . 
Inset: Replot of data collapse in terms of the variable LT Vph . 
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FIG. 6: Variation of the fitting error a against T c , for the 
total vortex correlation Xv according to Eq. [6l with B = 0.04. 
Inset: log-log fit assuming T c — 0. 

collaps o 29 ' 30 . The standard finite-size scaling expansion 
near T c , 

G(x) = a a + a\x + a 2 x 2 + . . . , (8) 

truncated to low order (0(x 5 )), is used as a smooth inter- 
polation function, which is convenient both for numerical 
data fitting and for providing a measure of the data col- 
lapse as the error denned by Eq. ([7]). Then we determine 
v from the best least-square fit of £(L, T)/L to this func- 
tion with x = TL X I V . Data for the largest temperatures 
are successively removed from the data collapse if this 
leads to a decrease of the fitting error, since in this case 
such data are presumably outside the asymptotic scaling 
regime where Eq. ([5]) applies. The best data collapse 
is obtained for temperatures T < 0.16. Fig. 7 shows 
that indeed the data for the phase variables £ p /, satisfy 
the scaling form with an exponent v v h — 1.29(15). Sim- 
ilarly, for the vortex variables, we obtain from the dada 
collapse in Fig. 8 the critical exponent v v — 1.03(25), 



although the data collapse is not as good. These expo- 
nents agree with each other within the estimated errors. 
It should be noted, however, that the scaling behavior in 
Eq. (O has only been verified here in a limited range of 
temperatures since there is no available data below 7/, 
due to the lack of equilibration. This can be seen more 
clearly by reploting the data collapse as L/£(L,T) ver- 
sus y = LT V , as presented in the insets of Figs. 7 and 
8. Indeed, L/£(L,T) tends to a linear behavior for large 
y as expected but the finite limit for y — > can not be 
explicitly verified. 

The zero-temperature transition, which was inferred 
here from the behavior of correlations at finite temper- 
atures as discussed above, is also consistent with the 
ground-state properties of the model. As first pointed 
out by Shih and Stroud^, the ground state is highly de- 
generate with an infinite number of different vortex con- 
figurations. Therefore, correlations at zero temperature, 
obtained by averaging over the different configurations, 
can decay to zero for large distances. One possibility is 
that such decay obeys a power law. The triangular an- 
tiferromagnetic Ising model^i, for example, display such 
behavior. Another possibility is that correlations decay 
exponentially as in the antiferromagnetic Ising model on 
a Kagome lattice 35 . The T c — scenario for the present 
model, favors a power-law decay of vortex correlations. 
In fact, the correlation function exponent r/ v associated 
with the transition can be obtained from the scaling law 
7 = v(2 — r\) and the above estimates of j v and v v , giving 
tj v = 0.90, which is greater than zero implying a decay 
of the correlation as C v (r) oc r^'' 11 , at zero temperature. 
To verify this behavior more quantitatively, we have also 
performed additional calculations of the finite-size de- 
pendence of the total correlation functions \ph and Xv 
near T = 0. However, because full equilibration was only 
possible for temperatures above the dynamical freezing 
transition Tf with the present Monte Carlo method, we 
have to rely on the approximation of using low-energy 
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FIG. 8: Scaling plot according to Eq. (J5|, assuming T c = 0, 
for the vortex correlation length £„, giving the estimate v v — 
1.03(25). Dashed line is a fit to the data with Eq. (JHJ) . Inset: 
Replot of the data collapse in terms of the variable LT Vv . 
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FIG. 9: Finite-size behavior of the total correlation function 
of phase Xph and vortex \ v variables obtained from low-energy 
states, using simulated annealing. Lines correspond to log-log 
fits, which give the estimates rjp^ — 0.4(2) and rj v = 0.4(3) 



states close to the ground state to perform the config- 
uration averages. The low-energy states were obtained 
using the simulated-annealing method 28 by gradual an- 
nealing from initial temperatures T > Tf. Fig. 9 shows 
a log-log plot of the finite -size behavior of X.ph 

and Xv 

obtained by averaging over 800 low-energy sates within 
a range of 0.16% above the known ground-state energy^ 
Eg = —1.2071. If correlations decay as a power-law in the 
ground state, then the finite-size behavior \ oc L 2 ~ n is 
expected near T = 0, for sufficiently large systems. This 
behavior is consistent with the data in Fig. 9 and from 
a log-log fit we estimate rj p h — 0.4(2) and r] v = 0.4(3), 
which are compatible within the estimated errors with 
the corresponding values obtained from 7 and u, using 
the finite-temperature scaling analysis. 

A particular feature of the frustrated XY model on 
the honeycomb lattice is the structure of the dual lattice, 
where the chiralities are defined, which has the form of 
a triangular lattice. Since the chiralities are Ising-like 
variables with antiferromagnetic interactions, if interac- 
tions further than nearest-neighbors are neglected, there 
is a geometric frustration effect which may appear simi- 
lar to the case of the triangular antiferromagnetic Ising 
model. The exact solution of this model^i shows that 
there is no phase transition at nonzero temperatures ( 
T c = 0) , which is the same scenario we find in the present 
case. However, while the correlation length for this Ising 
model diverges exponentially^, £ oc e 2 / T , corresponding 
to v — > 00 in Eq. ([5]), in the present case we find v v ~ 1, 
suggesting a power-law correlation and correspondingly 
different ground state properties. 

The near agreement of our estimates of the critical ex- 
ponents for vortex variables v v and phase variables u p h, 
assuming power-law divergent correlations, suggests that 
the T c = critical behavior may be described by a single 
divergent length scale. In a zero-temperature transition 
scenario where there is no decoupling, phase and vortex 
variables remain coupled on large length scales and one 



expects that the corresponding correlation lengths should 
diverge with a common critical exponent v p h = v v = v. It 
is interesting to note that this behavior is in contrast with 
that found for the two-dimensional XY model at irra- 
tional frustration on a square lattice^, where a decoupled 
zero-temperature transitional was found with significant 
different exponents v v >> v p h- 

An additional important feature of the zero- 
temperature transition found from the above analysis, 
is the temperature dependence of the relaxation time for 
phase and vortex equilibrium fluctuations. If T c — then 
the divergence of the relaxation time r as T — > T c is de- 
termined by thermal activation*!. Thus, we expect that 
t should increase exponentially with the inverse of tem- 
perature, corresponding to a dynamical exponent z — > 00 
in the usual power-law assumption r oc \T — T c \~ zv . To 
verify this behavior, we have in addition calculated the 
relaxation time r for different temperatures from the au- 
tocorrelation function of phase and vortex fluctuations, 
C p h{t) and C v (t), as 

1 f°° 

T ph,v = y, / dt C p h,v(t). (9) 

In these calculations, the starting configurations were 
taken from the equilibrium configurations obtained with 
the parallel tempering method and the subsequent time 
dependence was obtained from standard MC simulations 
at each fixed temperature. The results shown on the 
log-linear plot in Fig. 10 are indeed consistent with an 
activated behavior of r p h and t v . The straight lines in the 
plot indicate that the data can be fitted to an Arrhenius 
behavior, with temperature independent energy barriers 
E p h = 1-00 and E v = 1.02. In general, the energy barrier 
can be temperature dependent, scaling with the correla- 
tion length^! as E oc . The observed Arrhenius behav- 
ior indicates that these additional exponents are ip p h ~ 0, 
ip v « 0. The behavior for the relaxation time t p ^ is 
particularly important when the frustrated XY model 
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FIG. 10: Temperature dependence of the relaxation times of 
phase fluctuations r p h and vortex fluctuations r„ for system 
size L = 32. 

is applied to superconductors since it determines the 
temperature dependence of the linear resistivity ^ 31 ' 32 , 
Pl on 1/rph- As a consequence, pL should be finite at 
any nonzero temperature but decrease exponentially as 
temperature vanishes. 

V. CONCLUSIONS 

We have investigated the critical behavior of the fully 
frustrated XY model on a two-dimensional honeycomb 
lattice by Monte Carlo methods and a scaling analy- 
sis of the phase and vortex correlations. No evidence 
of an equilibrium phase transition is found at nonzero 
temperatures, in agreement with the conclusion of ref. 
I20I for a similar model. The absence of vortex order- 
ing at finite temperatures is also in agreement with the 
estimates of ref. |2l], which find that this can only be ob- 
served for very large system systems (L > 10 5 ). However, 
our finite-size scaling analysis is consistent with a zero- 
temperature transition, where the critical temperature 
vanishes and phase and vortex correlation lengths, t; p h 
and £„, diverge for decreasing temperatures as a power 
law with a common critical exponent, suggesting a cou- 
pled T c = transition scenario. Since both correlation 



lengths remain finite in the temperature range T > 0.11, 
where a KT transition^ or a spin-glass transition^ was 
proposed to take place in earlier MC simulations, our re- 
sults also indicate that these apparent transitions should 
be attributed to slow dynamics effects and not to an 
equilibrium phase transition. Whether a vortex ordering 
transition can, nevertheless, occur for very large system 
sizes (L > 10 5 ), as predicted in Ref. [2l|, can not be tested 
by our numerical simulation, which needs small system 
sizes in order to insure full equilibration. The main fea- 
tures of the correlation-length behavior obtained for the 
present model are the same as those found in other frus- 
trated XY models where the phase- coherence tempera- 
ture vanishes, such as, the two-dimensional gauge glass 
model 3 — and the XY model with irrational frustration 7 , 
but with different critical exponents. When applied to 
superconductors, the divergent correlation length in 
these models determine both the linear an nonlinear re- 
sistivity behavior, leading to a thermally activated lin- 
ear resistivity and nonlinear current-voltage scaling at 
low temperatures ^ 31 ' 32 . Similar behavior should be ob- 
served for Josephson-junction arrays on a honeycomb lat- 
tice and superconducting films with a triangular pattern 
of nanohole o 16 ' 39 in a perpendicular magnetic field corre- 
sponding to half flux quantum per plaquette, since both 
systems can be modeled by a frustrated XY model on a 
honeycomb lattice. The predicted thermal activated be- 
havior for the linear resistivity seems to be already been 
observed in the latter system^. 
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